class OrderedProbitNaive(pm.Categorical):
    def __init__(self, eta, cutpoints, *args, **kwargs):
        self.eta = tt.as_tensor_variable(pm.floatX(eta))
        self.cutpoints = tt.as_tensor_variable(cutpoints)
        pa = pm.distributions.dist_math.std_cdf(self.cutpoints - tt.shape_padright(self.eta))
        p_cum = tt.concatenate(
            [
                tt.zeros_like(tt.shape_padright(pa[..., 0])),
                pa,
                tt.ones_like(tt.shape_padright(pa[..., 0])),
            ],
            axis=-1,
        )
        p = p_cum[..., 1:] - p_cum[..., :-1]
        super().__init__(p=p, *args, **kwargs)

n1_c = 300; n2_c = 300; n3_c = 300
cluster1 = np.random.randn(n1_c) + -1
cluster2 = np.random.randn(n2_c) + 0
cluster3 = np.random.randn(n3_c) + 2

x = np.concatenate((cluster1, cluster2, cluster3))
y = np.concatenate((1*np.ones(n1_c),
                    2*np.ones(n2_c),
                    3*np.ones(n3_c))) - 1

# Add an outlier point
x[0] = -10
y[0] = 2

# Ordered logistic regression
with pm.Model() as model:
    cutpoints = pm.Normal("cutpoints", mu=[-1,1], sigma=10, shape=2,
                          transform=pm.distributions.transforms.ordered)
    y_ = OrderedProbitNaive("y", cutpoints=cutpoints, eta=x, observed=y)
    tr = pm.sample(1000, chains=1)